if (class(model) %in% "character") {
  model <- ModelTBBCGEngland::read_libbi(file.path(model_dir, "libbi", "posterior"))
}

priors <- readRDS(file.path(model_dir, "data", "prior-params.rds"))
model
## Wrapper around LibBi
## ======================
## Model:  Baseline 
## Run time:  406.828  seconds
## Number of samples:  1000 
## State trajectories recorded:  E H L P S T_E T_P YearlyAgeCases YearlyCases YearlyEAgeCases YearlyECases YearlyEPulCases YearlyEPulDeaths YearlyNonUKborn YearlyPAgeCases YearlyPCases YearlyPulCases YearlyPulDeaths 
## Observation trajectories recorded:  YearlyAgeInc YearlyHistPInc YearlyInc 
## Parameters recorded:  alpha_t_decay alpha_t_init c_eff c_hist chi_init delta epsilon_h_0_4 epsilon_h_15_89 epsilon_h_5_14 epsilon_l_0_4 epsilon_l_15_89 epsilon_l_5_14 HistMeasError kappa_0_4 kappa_15_89 kappa_5_14 M mu_e_0_14 mu_e_15_59 mu_e_60_89 mu_p_0_14 mu_p_15_59 mu_p_60_89 nu_e_0_14 nu_e_15_89 nu_p_0_14 nu_p_15_89 phi_0_14 phi_15_59 phi_60_89 rho_0_14 rho_15_59 rho_60_89 Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89 zeta_p_0_14 zeta_p_15_59 zeta_p_60_89
traces <- coda::mcmc(rbi::get_traces(model))

coda::rejectionRate(traces)
##               M           c_eff          c_hist           delta 
##        0.995996        0.995996        0.995996        0.995996 
##   epsilon_h_0_4  epsilon_h_5_14 epsilon_h_15_89       kappa_0_4 
##        0.995996        0.995996        0.995996        0.995996 
##      kappa_5_14     kappa_15_89   epsilon_l_0_4  epsilon_l_5_14 
##        0.995996        0.995996        0.995996        0.995996 
## epsilon_l_15_89        phi_0_14       phi_15_59       phi_60_89 
##        0.995996        0.995996        0.995996        0.995996 
##    Upsilon_0_14   Upsilon_15_59   Upsilon_60_89        rho_0_14 
##        0.995996        0.995996        0.995996        0.995996 
##       rho_15_59       rho_60_89       nu_p_0_14      nu_p_15_89 
##        0.995996        0.995996        0.995996        0.995996 
##       nu_e_0_14      nu_e_15_89     zeta_p_0_14    zeta_p_15_59 
##        0.995996        0.995996        0.995996        0.995996 
##    zeta_p_60_89     zeta_e_0_14    zeta_e_15_59    zeta_e_60_89 
##        0.995996        0.995996        0.995996        0.995996 
##       mu_p_0_14      mu_p_15_59      mu_p_60_89       mu_e_0_14 
##        0.995996        0.995996        0.995996        0.995996 
##      mu_e_15_59      mu_e_60_89        chi_init    alpha_t_init 
##        0.995996        0.995996        0.995996        0.000000 
##   alpha_t_decay   HistMeasError 
##        0.995996        0.995996
coda::autocorr.plot(traces)

- Evaluate Posterior Traces

plot_param(model, scales = "free", plot_type = "trace", burn_in = 0) + theme(legend.position = "none")

plot_param(model, prior_params = priors, scales = "free")

param_sum <- plot_param(model, prior_params = priors, plot_data = FALSE) %>% 
  group_by(Distribution, parameter, length) %>% 
  summarise(median = median(value), 
            lll = quantile(value, 0.025), 
            hhh = quantile(value, 0.975)) %>% 
  mutate(value = pretty_ci(median, lll, hhh, sep = ", ")) %>%
  group_by(Distribution, parameter) %>% 
  mutate(Parameter = case_when(max(length) > 1 ~ paste(parameter, 1:n(), sep = "-"),
                               TRUE ~ as.character(parameter))
         ) %>% 
  ungroup %>% 
  select(Distribution, Parameter, value) %>% 
  spread(key = "Distribution", value = "value") %>% 
  select(Parameter, Prior, Posterior)

saveRDS(param_sum, file.path(model_dir, "data", "sum_prior_posterior.rds"))
knitr::kable(param_sum)
Parameter Prior Posterior
alpha_t_decay -0.13 (-0.24, -0.03) -0.19 (-0.26, -0.05)
alpha_t_init 0.84 (0.76, 0.90) 0.50 (0.11, 0.89)
c_eff 2.46 (0.14, 4.88) 3.84 (1.43, 4.18)
c_hist 12.58 (10.12, 14.89) 12.00 (10.81, 13.45)
chi_init 0.19 (0.08, 0.29) 0.23 (0.18, 0.30)
delta 0.78 (0.70, 0.85) 0.78 (0.77, 0.84)
epsilon_h_0_4 1.53 (0.10, 4.88) 1.11 (0.69, 3.72)
epsilon_h_15_89 25.72 (1.45, 77.21) 9.04 (4.43, 70.81)
epsilon_h_5_14 3.72 (0.16, 11.99) 5.85 (4.21, 7.45)
epsilon_l_0_4 594.41 (29.25, 1714.03) 325.02 (183.07, 1257.01)
epsilon_l_15_89 1066.55 (38.46, 3368.79) 44.09 (31.10, 51.12)
epsilon_l_5_14 537.07 (31.69, 1541.71) 837.75 (28.05, 1249.42)
HistMeasError 1.00 (0.60, 1.37) 1.16 (1.03, 1.38)
kappa_0_4 0.88 (0.05, 2.78) 1.01 (0.06, 1.21)
kappa_15_89 1.11 (0.04, 3.47) 1.26 (0.62, 1.87)
kappa_5_14 0.97 (0.05, 3.10) 0.61 (0.18, 0.73)
M 0.25 (0.01, 0.49) 0.21 (0.01, 0.40)
mu_e_0_14 275.27 (211.92, 336.54) 256.76 (249.02, 301.04)
mu_e_15_59 194.13 (69.41, 321.66) 160.81 (127.82, 282.82)
mu_e_60_89 36.50 (1.43, 106.56) 35.09 (0.80, 42.95)
mu_p_0_14 240.25 (154.15, 330.24) 248.58 (195.58, 270.37)
mu_p_15_59 88.24 (5.27, 259.49) 41.89 (14.04, 68.54)
mu_p_60_89 46.99 (2.35, 148.94) 75.42 (41.67, 127.12)
nu_e_0_14 0.17 (0.00, 1.30) 0.15 (0.04, 0.26)
nu_e_15_89 0.33 (0.01, 1.67) 0.71 (0.08, 2.72)
nu_p_0_14 0.13 (0.00, 0.71) 0.07 (0.01, 0.36)
nu_p_15_89 0.23 (0.01, 1.17) 0.16 (0.06, 0.52)
phi_0_14 0.58 (0.30, 1.06) 0.47 (0.40, 0.55)
phi_15_59 0.62 (0.30, 1.15) 0.38 (0.20, 0.78)
phi_60_89 0.60 (0.29, 1.12) 0.68 (0.31, 0.97)
rho_0_14 0.30 (0.26, 0.34) 0.29 (0.27, 0.33)
rho_15_59 0.65 (0.64, 0.66) 0.65 (0.64, 0.65)
rho_60_89 0.54 (0.52, 0.55) 0.53 (0.53, 0.54)
Upsilon_0_14 0.63 (0.61, 0.65) 0.64 (0.62, 0.64)
Upsilon_15_59 0.71 (0.70, 0.71) 0.71 (0.70, 0.71)
Upsilon_60_89 0.75 (0.74, 0.76) 0.75 (0.75, 0.76)
zeta_e_0_14 123.74 (58.77, 187.11) 174.99 (106.52, 193.46)
zeta_e_15_59 62.71 (3.84, 175.77) 88.10 (5.55, 156.23)
zeta_e_60_89 132.08 (55.59, 209.75) 90.31 (84.13, 191.14)
zeta_p_0_14 97.26 (17.74, 184.18) 89.02 (51.59, 107.26)
zeta_p_15_59 81.34 (3.53, 252.31) 133.28 (27.51, 146.85)
zeta_p_60_89 122.28 (16.92, 246.71) 155.70 (10.61, 231.57) 
plot_state(model, summarise = TRUE)

plot_state(model, summarise = TRUE, start_time = 59)

plot_state(model, states = c("YearlyHistPInc", "YearlyInc"), summarise = TRUE, start_time = 49) -> p2

p2

pred_obs <- plot_state(model, states = c("YearlyHistPInc", "YearlyInc", "YearlyAgeInc"), summarise = FALSE, start_time = 59, plot_data = FALSE)

obs <- ModelTBBCGEngland::setup_model_obs() %>% 
  bind_rows(.id = "state") %>% 
  mutate(time = time + 1931)

pred_obs <- pred_obs %>% 
  dplyr::filter(Average == "median") %>% 
  mutate(pred_value = pretty_ci(Count, lll, hhh, digits = 0)) %>% 
  left_join(obs, by = c("state" = "state", "time" = "time", "age" = "age")) %>% 
  select(Observation = state, Age = age, Year = time,  `Observed Incidence` = value, `Predicted Incidence` = pred_value) %>% 
  mutate(Year = Year)


sum_obs_table <- pred_obs %>% 
  dplyr::filter(Observation %in% c("YearlyInc", "YearlyHistPInc")) %>% 
  drop_na(`Observed Incidence`) %>% 
  select(-Age) %>% 
  mutate(Observation = case_when(Observation == "YearlyInc" ~ "All TB cases",
                                 Observation == "YearlyHistPInc" ~ "Pulmonary TB cases"))


saveRDS(sum_obs_table, file.path(model_dir, "data", "sum_obs_table.rds"))

age_cases_table <- pred_obs %>% 
  dplyr::filter(Observation %in% c("YearlyAgeInc")) %>% 
  drop_na(`Observed Incidence`) %>% 
  select(-Observation) %>% 
  group_by(Year) %>% 
  mutate(Age = c(paste0(seq(0, 65, 5), "-", seq(4, 69, 5)), "70-89")) %>% 
  ungroup


saveRDS(age_cases_table, file.path(model_dir, "data", "age_cases_table.rds"))
knitr::kable(sum_obs_table)
Observation Year Observed Incidence Predicted Incidence
Pulmonary TB cases 1990 3618 4181 (2798 to 4349)
Pulmonary TB cases 1991 3596 4100 (2745 to 4400)
Pulmonary TB cases 1992 3816 3905 (2803 to 4318)
Pulmonary TB cases 1993 3961 3790 (2604 to 4275)
Pulmonary TB cases 1994 3694 3731 (2558 to 4103)
Pulmonary TB cases 1995 3711 3709 (2614 to 4067)
Pulmonary TB cases 1996 3690 3437 (2451 to 4033)
Pulmonary TB cases 1997 3947 3503 (2431 to 3988)
Pulmonary TB cases 1998 3926 3304 (2434 to 3910)
Pulmonary TB cases 1999 3827 3286 (2384 to 3807)
All TB cases 2000 1803 2417 (1618 to 3261)
All TB cases 2001 1866 2319 (1521 to 3107)
All TB cases 2002 1833 2105 (1521 to 3102)
All TB cases 2003 1685 2124 (1371 to 3090)
All TB cases 2004 1776 1949 (1290 to 3058) 
knitr::kable(age_cases_table)
Age Year Observed Incidence Predicted Incidence
0-4 2000 75 1 (0 to 14)
5-9 2000 59 0 (0 to 12)
10-14 2000 75 1 (0 to 21)
15-19 2000 112 4 (2 to 35)
20-24 2000 133 6 (2 to 46)
25-29 2000 116 12 (6 to 89)
30-34 2000 147 28 (21 to 108)
35-39 2000 99 61 (48 to 141)
40-44 2000 83 117 (80 to 234)
45-49 2000 105 215 (123 to 318)
50-54 2000 116 296 (175 to 369)
55-59 2000 112 344 (258 to 500)
60-64 2000 101 440 (316 to 457)
65-69 2000 99 462 (321 to 498)
70-89 2000 371 374 (260 to 381)
0-4 2001 113 2 (0 to 23)
5-9 2001 65 1 (0 to 13)
10-14 2001 51 1 (0 to 20)
15-19 2001 130 1 (1 to 46)
20-24 2001 145 5 (2 to 64)
25-29 2001 127 9 (8 to 80)
30-34 2001 122 29 (15 to 112)
35-39 2001 128 45 (37 to 169)
40-44 2001 107 116 (62 to 220)
45-49 2001 101 170 (113 to 306)
50-54 2001 96 281 (165 to 386)
55-59 2001 91 331 (233 to 469)
60-64 2001 94 406 (253 to 482)
65-69 2001 105 472 (298 to 476)
70-89 2001 391 364 (280 to 395)
0-4 2002 102 1 (0 to 16)
5-9 2002 62 0 (0 to 19)
10-14 2002 64 2 (0 to 15)
15-19 2002 129 2 (0 to 45)
20-24 2002 148 8 (3 to 59)
25-29 2002 120 16 (7 to 87)
30-34 2002 112 27 (14 to 116)
35-39 2002 127 62 (33 to 178)
40-44 2002 98 102 (59 to 219)
45-49 2002 89 174 (113 to 307)
50-54 2002 104 241 (176 to 398)
55-59 2002 116 330 (203 to 468)
60-64 2002 88 411 (270 to 452)
65-69 2002 90 435 (315 to 483)
70-89 2002 384 345 (265 to 348)
0-4 2003 74 2 (0 to 20)
5-9 2003 46 1 (1 to 31)
10-14 2003 59 2 (0 to 27)
15-19 2003 102 4 (1 to 44)
20-24 2003 116 6 (4 to 59)
25-29 2003 137 11 (4 to 83)
30-34 2003 118 32 (10 to 123)
35-39 2003 113 53 (27 to 177)
40-44 2003 109 98 (47 to 228)
45-49 2003 83 132 (87 to 304)
50-54 2003 102 230 (126 to 359)
55-59 2003 92 358 (206 to 424)
60-64 2003 98 378 (259 to 456)
65-69 2003 94 418 (297 to 453)
70-89 2003 342 339 (241 to 379)
0-4 2004 112 3 (0 to 28)
5-9 2004 71 2 (0 to 25)
10-14 2004 81 2 (1 to 29)
15-19 2004 111 3 (1 to 54)
20-24 2004 120 9 (2 to 74)
25-29 2004 136 11 (5 to 99)
30-34 2004 129 22 (8 to 126)
35-39 2004 103 50 (23 to 154)
40-44 2004 128 93 (43 to 226)
45-49 2004 94 140 (82 to 240)
50-54 2004 103 208 (135 to 340)
55-59 2004 98 291 (185 to 429)
60-64 2004 86 371 (256 to 414)
65-69 2004 92 414 (311 to 434)
70-89 2004 312 338 (212 to 369) 
plot_state(model, states = c("YearlyAgeInc"), start_time = 59) + theme(legend.position = "none")